#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double *u, *temp, *r;
int n,steps,test_n;
double L,dx;
double J,J2,M;
double alpha, beta;

#define _get(x,i,j,k) (x)[((i)*n+(j))*n+(k)]
#define prev(i) (((i)==0)?1:(i)-1)
#define get(x) _get(x,i,j,k)

#define X (i*dx)
#define Y (j*dx)
#define Z (k*dx)

#define nx _get(u,i+1,j,k)
#define ny _get(u,i,j+1,k)
#define nz _get(u,i,j,k+1)
#define px _get(u,prev(i),j,k)
#define py _get(u,i,prev(j),k)
#define pz _get(u,i,j,prev(k))
#define ppx _get(u,i-1,j,k)
#define ppy _get(u,i,j-1,k)
#define ppz _get(u,i,j,k-1)
#define U get(u)

double pow8(double x){
    x = x*x;
    x = x*x;
    x = x*x;
    return x;
}

#define F (15.75*J2*(X*X+Y*Y)/pow8((U+1)*get(r)+M))

void update_u(){
    int i,j,k;
    for(i=0; i<n-1; i++)
	for(j=0; j<n-1; j++)
	    for(k=0; k<n-1;k++)
		get(temp) = alpha*(nx+ny+nz+px+py+pz)/6
		            +U*(1-alpha*(1+F*dx*dx/6));
    k = n-1;
    for(i=0; i<n-1; i++)
	for(j=0; j<n-1; j++)
	    get(temp) = (4*pz - ppz - 2*dx*U - X*(nx - px) - Y*(ny - py)) / 3;

    i = n-1;
    for(j=0; j<n-1; j++)
	for(k=0; k<n-1; k++)
	    get(temp) = (4*px - ppx - 2*dx*U - Y*(ny - py) - Z*(nz - pz)) / 3;
    
    j = n-1;
    for(k=0; k<n-1; k++)
	for(i=0; i<n-1; i++)
	    get(temp) = (4*py - ppy - 2*dx*U - Z*(nz - pz) - X*(nx - px)) / 3;

    i = j = n-1;
    for(k=0; k<n-1; k++)
	get(temp) = (4*(px+py) - ppx - ppy - 2*dx*U - Z*(nz - pz)) / 6;

    j = k = n-1;
    for(i=0; i<n-1; i++)
	get(temp) = (4*(py+pz) - ppy - ppz - 2*dx*U - Z*(nx - px)) / 6;

    k = i = n-1;
    for(j=0; j<n-1; j++)
	get(temp) = (4*(pz+px) - ppz - ppx - 2*dx*U - Z*(ny - py)) / 6;

    i = j = k = n-1;
    get(temp) = (4*(px+py+pz) - ppx - ppy - ppz - 2*dx*U) / 9;

    for(i=0; i<n; i++)
	for(j=0; j<n; j++)
	    for(k=0; k<n;k++)
		U = get(temp);
}

void get_u(){
    J2 = J*J;
    n = L / dx;
    printf("L=%lf, dx=%e, n=%d, J=%lf, M=%lf\n", L, dx, n, J, M);
    u = malloc(sizeof(double)*n*n*n);
    r = malloc(sizeof(double)*n*n*n);
    for(int i=0; i<n; i++)
	for(int j=0; j<n; j++)
	    for(int k=0; k<n; k++){
		get(r) = sqrt(X*X+Y*Y+Z*Z);
		U = exp(-X*X-Y*Y-Z*Z);
		//U = 1;
	    }
    temp = malloc(sizeof(double)*n*n*n);
    for(int step=0; step<steps; step++){
	update_u();
	if(step%test_n==0){
	    double A=0, B=0;
	    double Ax,Ay,Az,Bx,By,Bz;
	    for(int i=0;i<n;i++)
		for(int j=0; j<n; j++)
		    for(int k=0; k<n; k++){
			if(U>A){
			    A  = U;
			    Ax = X;
			    Ay = Y;
			    Az = Z;
			}
			if(U<B){
			    B  = U;
			    Bx = X;
			    By = Y;
			    Bz = Z;
			}
		    }
	    printf("step:%d max:%e at (%lf, %lf, %lf), min:%e at (%lf, %lf, %lf)\n",
		    step,A,Ax,Ay,Az,B,Bx,By,Bz);
	}
    }
}

